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1 . Formulation of a Fast Iterative Method 
f o r Transonic Flow Calculations 


Reliable but slow methods for calculating transonic 
flows have been developed in recent years [1,2,3]. These » 

use central difference formulas in the subsonic zone and 
upwind difference formulas in the supersonic zone to ensure 
the proper region of dependence and jump conditions. The 
resulting difference equations are then solved by an itera- 
tion procedure derived from the method of successive 
overrelaxation. This note describes results obtained by 
using a fast Poisson solver to accelerate the rate of 
convergence of the iterative scheme. 


Note: This work was also partially supported by NASA 

Grants NGR-33016-167 and -201. 
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It was proposed by Martin and Lomax (4] that a fast 
elliptic solver could be used to generate an iterative 
scheme for solving the difference equations appearin j in 
compressible flow calculations. In the simplest case 
consider the small disturbance equation 

2 

( 1-M )<J> + <P - 0 

,v xx T yy 

where $ is the velocity potential, and M is the local 
Mach number, which is related to the free stream Mach 
number M by the formula 

M 2 - Mf(l + (Y+l)$ x > • 

An iterative scheme can be constructed by putting the 
Laplacian on the left and the nonlinear terms on the right. 
Let v^ be the solution for <(> at the nth iteration. Then 


Av 


n+1 


K 2 X V 


3 x 


n 


To see that this can be expected to converge consider the 

2 2 

linearized equation where M is replaced by M^. If P 
and Q are nonnegative finite difference operators represent- 


ing - — y and 
3x 


_3 

3y 


we have 


(P+Q) v n+ i - M Pv n 


or 


_l/2 
P ' v 


n+1 


.,2 v 1/2 
M K P ' v_ 


n 


where 


K = P 1/2 (P+Q) _1 P ' ' 
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Th us 


lp 1/2 v < M 2 lKl I r 1/2 V I 

n+1 — n 


and since K is Hermitian 


I Kl = A (K) 
max 


(x,Kx) 
X (x, x) 


(y ,Pv ) 

max fy , Pyf" +’Ty”7Qy) 


where 


n l /2 „ 1/2 

P ' y = K x 


Thus 


«P 1/2 v , , I < M 2 IP 1/2 v I 


n+1 - 


n 


This estimate serves to indicate that for subsonic flows 
the scheme should converge at a rate independent of the 
mesh size. 

The above analysis also suggests what it is doubtful 
whether such a scheme would converge for a supersonic flew 
with M > 1. The argument presented in the Appendix 
in fact indicates that the scheme would uefinitely not 
converge for a linearized supers oi.ic flow. It thus 
appears that the fast elliptic solver used on its own is 
not likely to lead to good conver ence when the supersonic 
zone is large. If, however, it could be supplemented with 
another scheme which give fast convergence in the super- 
sonic zone, the iwo in combination might produce an 
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effective iterative scheme. In fact the scandard line 
relaxation method for transonic flow calculations is such 
a scheme. For the small disturbance equation, or in the 
case cf the full potential equation with the flow aligned with 
one coordinate direction, the method consists of freez- 

ing the nonlinear coefficients at values determined from 
the previous iteration and solving the resulting wava 
equation in the supersonic zone by a marching procedure. 

Thus an exact solution of the supersonic zone could be 
obtained in one step if the correct coefficients and data 
at the sonic line could be inserted. 

Thus the following scheme is proposed: use a two 

stage iteration, in which the first stage is a step using 
the Laplacian on the left-hand side, and the second stag*! 
consists of a fixed number of relaxation steps to stabilize 
the supersonic zone. For the case of flow aligned with 
the coordinate system a single relaxation step should be 
sufficient. The full potential equation in a curvilinear 
coordinate systom with an arbitrary flow direction requires 
the use of a rotated upwind difference scheme [3]. Then a 
simple marching scheme can no longer be used in the super- 
sonic zone, and .everal relaxation steps may be required 
in the second stage. 


2 . Application to the Transonic Potential Flow Equation in 






a Mapped Domain 

The use of a fast Poisson solver requires a simple domain 
such as a rectangle. This leads to a difficulty in applying 
the proposed method to an exterior flow problem with an infinite 
domain. This can be circumvented by using the full potential 
equation and mapping the exterior of the profile onto the 
interior of a circle. If 2 ttE is the circulation it is conveni- 
ent to use a reduced potential G defined by 

4> «* G + * - - E 6. 

Then G is finite and single valued. Now a fast solver for 
Poisson's equation in polar coordinates can be used in the 
first stage of the iteration. For this purpose a scheme using 
the Buneman algorithm in the 0 direction has been programmed. 

Two variants of this approach have been tried. The first 
treats the potential equation in quasilinear form. The residual 
at each point is evaluated as 

R * (a 2 -v 2 ) r (r G r ) + (a 2 -u 2 ) G, )0 - 2uv r(G r0 +G 0 - E) 

+ (u 2 -v 2 ) r G r + (u 2 +v 2 ) H (] + v H r ) , 

where H is the modulus of the transformation to the exterior 

of the circle, u and v are the velocity components 

2 

r(G fl -E) - si.i 0 r G - cos 0 

u .. r 


and a is the speed of sound. If y 


is the ratio of 


3 po c i f i c heats, a is determined from the stagnation speed 


of sound a Q by the relation 


2 

a 


2 

a 0 



2 2 
(u%v ) 


In evaluating R upwind differencing is used in the usual 
manner at supersonic points. In the first stage of the 
iteration the correction C is determined by solving 

r h (r V * c ee - 

and then G is updated by the rule 

G + ■ G + to C 


where the superscript + denotes the new value, and 
a) is an overrelaxation factor. In the second stage of 
the iteration an ordinary relaxation step is used. 

Results with this approach have been quite promising. 
Numerical tests have confirmed that the scheme sometimes 
diverges when the relaxation step is not included. When it is 
included fast convergence has been obtained. Figures 1 
and 2 show typical results. In each calculation the 
calculation was performed first on a mesh with 64 cells in 
the 0 direction and 16 cells in the r direction, and then 
on a mesh with 128 * 32 cells. The interpolated coarse mesh 
solution was used as the starting point for the fine mesh 
calculation. The largest absolute value of the residual 
anywhere in the field was used as a measure of convergence. 
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The first example is the 64A410 airfoil at Mach .720. 

>1 -9 

In this case the residual was reduced from ''•lO to 10 

- 3 .9 

in 26 cycles on the coarse mesh, and then from ^10 to 10 

in 21 cycles on the fine mesh, each cycle consisting of 

one Poisson step plus one relaxation \tep. The Poisson 

step takes about the same time a*' 2 relaxation steps, so 

each complete cycle requires about the same time as 3 

relaxation steps. On the CDC 6600 at the ERDA Computing 

Facility at New York University one complete cycle on the 

fine mesh takes about 1.5 seconds. The entire calculation 

for the 64A410 took 48 seconds. The second example shows 

a shock-free supercritical airfoil designed by Garabedian (61 . 

In this case 27 cycles were required to reduce the largest 
-9 

residual to 10 on the coarse mesh, and another 27 cycles 

-9 

to reduce it to 10 on the fine mesh. In corresponding 
calculations using relaxation steps without the Poisson 
steps the largest residual was still v-lO 6 after ’000 cycles 
on the fine mesh. 

The second variant treats the potential equation in 
conservation form using a rotated difference scheme in the 
supersonic zone [5]. In this case the residual is evaluated 
as 

K = r f r - (rpV) + <pu> 


where 


» ■ G e * E - ^ ■ 


V = r G r - , 


and 


is the density. If M is the free stream Mach 
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number p is determined from the speed of sound a 
by the relation 


.Y-i 


M 


2 2 


Wow in the first stage of the iteration the correction C 
is determined by solving 


fr <r 


C r» 


+ c 


60 


R 

P 


The second stage consists of k relaxation cycles. Typically 
k » .. . 

In this case also the combined iteration has proved 
to give faster convergence than the simple relaxation 
method. The improvement is not as great as with the 

first variant, however, because of the need to use more 
relaxation steps than Poisson steps. As an example of 
the application of the method, Figure 3 shows the pressure 
distribution for the 64A410 at Mach .720 recalculated using 
conservation form. The proper theoretical jump condition 
is now satisfied, as can be seen. In this calculation 
the number of cycles required to reduce the largest residual 
to 10 wac 37 on the coarse mesh and 40 on the fine mesh. 
Each cycle consisted of 1 Poisson step plus 5 relaxation 
steps, and took about the same amount of time as 7 relaxa- 
tion steps, so the fine mesh calculation is equivalent to a 
little under 300 relaxation steps, which would be enough 
to reduce the largest residual to 'vio” 5 using relaxation 
alone. 
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3 . 


Conclusion 


It is concluded that the use of a fast elliptic solver 
in combination with relaxation is an effective way to 
accelerate the convergence of transonic flow calculations, 
particularly when a marching scheme can be used to treat 
the supersonic zone in the relaxation process. Other 
methods of preventing divergence in the supersonic zone 
should be investigated. Possibly it would be sufficient 
to sweep only the supersonic zone in the second stage of 
the iteration. This would lead to further useful savings 
of computer time. 
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NACA 64A410 AIRFOIL 

MACH MO. .720 ANGLE OF ATTACK 0° 

LIFT COEFFICIENT .6289 DRAG COEFFICIENT .0040 

Figure 1 
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78-06 AIRFOIL 


MACH NO. .780 ANGLE OF ATTACK 0° 

LIFT COEFFICIENT .5867 DRAG COEFFICIENT .0006 

Figure 2 
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- 1.20 - 1.60 
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NACA 64A410 AIRFOIL 

MACH NO. .720 ANGLE OF ATTACK 0° 

LIFT LIFT COEFFICIENT .6654 DRAG COEFFICIENT .0029 


Figure 3 
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Appendix. Analysis of the Poisson Iteration for the 


Linearized Equation with M > 1. 


Let the equation 


(1-M 2 )4> - 0 , 

' r xx yy ' 

, 2 

with M a constant > 1, be approximated with equal mesh 
spacing in the x and y directions by the Murman difference 
scheme 


(l-M 


i#j+l *ij 


in which an upwind difference formula is used for <J> . 

X X 

Denoting updated values by the superscript +, consider the 
iteration 


<b , . , “ 2d,, + a . , , + d> . . , — 2d 1 .. + <P . , , 

^i+l # j x} i“l#j i»D+l ID i»j-l 


i+1 


— 2 <p , . + <p , , 

,j ij i-l » 


+ (M _1) ( *ij' 2 *i-l,j + *i-2,j ) - 


Let i and j both run from 1 to n and define the n x n matrix 


•2 1 
1 -2 

1 


1 

-2 1 
1-2 1 


Also define the n x n matrix 
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R » 


1 

-2 1 
1-2 1 
1-2 1 


and let 4> be the macrix with entries 
iteration can be written as 

T 4> + + <t> + T ■ (T + aR) 4> 


whe re 


a 


M 2 -l > 


0 . 


Then the 


Also it is easily verified that 

R = T S 


where 



-l/(n-H) 
-2/ (n+1) 
- 3/ (n+1) 
-4/ (n+1) 


Thus 


T <t> + + 4> + T = T (I + aS) ♦ 


If 


X <J> 


then <I> (with its elements suitably ordered to form a vector) 
is an eigenvector of the iteration matrix. Consider the 
form 

T 

$> = uv 
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where v is an eigenvector of T 


T v * p v 


Then $ is an eigenvector if 


T ( I + aS) uv 


T T 

A(Tuv + uv T) 


This is satisfied if 


A (T + pi ) uv J 


where 


Pp u * A u 


-1 


P ( =* (T+pI) T ( I+aS) . 


Thus the eigenvectors of the iteration matrix can be 
expressed as 


u v 


where v is an eigenvector of T with eigenvalue p and 

u is an eigenvector of P ( , and the eigenvalues of the 

iteration matrix are the eigenvalues of P^ for 

P » p ^ » Pj » • • • / p ♦ For large n the smallest eigenvalue 

p of T is of order But as p ■* 0 the eigenvalues 

n 

of P approach those of 

P Q * I + aS 

and correspondingly some of the eigenvalues of the itera- 
tion matrix approach 

1 + a A. 


where A^ are the eigenvalues of S 


Now 


det (AI-S) 


. A n + n A n-1 
n+1 


♦ 


1 

n+l 


Thus if the polynomial 

J - (n+1) A n ♦ n A”” 1 ...-tl 

n 


has a root in the right half plane the corresponding 
eigenvalues of P Q will lie outside the unit circle. 
Applying the Routh Hurwitz test, it can be verified 
that J n has at least one root in the right half plane 
when n > 4 . 
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Government sponsored work. Neither the 
United States, nor the Administration, 
nor any person acting on behalf of the 
Admlnist nation : 

A. M akes any warranty or representation, 
express or Implied, with respect to the 
accuracy, completeness, or usefulness of 
the Information contained In this report, 
or that the use of any Information, 
apparatus, method, or process disclosed 
In this report may not Infringe privately 
owned rights; or 

B. Assumes any liabilities with respect to 
the use of, or for damages resulting from 
the use of any information, apparatus, 
method, or process disclosed in this 
report . 

As used in the above, "person acting on behalf 
of the Administration" Includes any employee 
or contractor of the Administration, or 
employee of such contractor, to the extent 
that such employee or contractor of the 
Administration, or employee of such contractor 
prepares, disseminates, or provides access to, 
any information pursuant to his employment or 
contract with the Administration, or his 
employment with such contractor. 
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